home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zsyr2k.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  10.7 KB  |  328 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        UPLO, TRANS
  8.       INTEGER            N, K, LDA, LDB, LDC
  9.       COMPLEX*16         ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  ZSYR2K  performs one of the symmetric rank 2k operations
  18. *
  19. *     C := alpha*A*B' + alpha*B*A' + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*A'*B + alpha*B'*A + beta*C,
  24. *
  25. *  where  alpha and beta  are scalars,  C is an  n by n symmetric matrix
  26. *  and  A and B  are  n by k  matrices  in the  first  case  and  k by n
  27. *  matrices in the second case.
  28. *
  29. *  Parameters
  30. *  ==========
  31. *
  32. *  UPLO   - CHARACTER*1.
  33. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  34. *           triangular  part  of the  array  C  is to be  referenced  as
  35. *           follows:
  36. *
  37. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  38. *                                  is to be referenced.
  39. *
  40. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  41. *                                  is to be referenced.
  42. *
  43. *           Unchanged on exit.
  44. *
  45. *  TRANS  - CHARACTER*1.
  46. *           On entry,  TRANS  specifies the operation to be performed as
  47. *           follows:
  48. *
  49. *              TRANS = 'N' or 'n'    C := alpha*A*B' + alpha*B*A' +
  50. *                                         beta*C.
  51. *
  52. *              TRANS = 'T' or 't'    C := alpha*A'*B + alpha*B'*A +
  53. *                                         beta*C.
  54. *
  55. *           Unchanged on exit.
  56. *
  57. *  N      - INTEGER.
  58. *           On entry,  N specifies the order of the matrix C.  N must be
  59. *           at least zero.
  60. *           Unchanged on exit.
  61. *
  62. *  K      - INTEGER.
  63. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  64. *           of  columns  of the  matrices  A and B,  and on  entry  with
  65. *           TRANS = 'T' or 't',  K  specifies  the number of rows of the
  66. *           matrices  A and B.  K must be at least zero.
  67. *           Unchanged on exit.
  68. *
  69. *  ALPHA  - COMPLEX*16      .
  70. *           On entry, ALPHA specifies the scalar alpha.
  71. *           Unchanged on exit.
  72. *
  73. *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  74. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  75. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  76. *           part of the array  A  must contain the matrix  A,  otherwise
  77. *           the leading  k by n  part of the array  A  must contain  the
  78. *           matrix A.
  79. *           Unchanged on exit.
  80. *
  81. *  LDA    - INTEGER.
  82. *           On entry, LDA specifies the first dimension of A as declared
  83. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  84. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  85. *           be at least  max( 1, k ).
  86. *           Unchanged on exit.
  87. *
  88. *  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
  89. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  90. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  91. *           part of the array  B  must contain the matrix  B,  otherwise
  92. *           the leading  k by n  part of the array  B  must contain  the
  93. *           matrix B.
  94. *           Unchanged on exit.
  95. *
  96. *  LDB    - INTEGER.
  97. *           On entry, LDB specifies the first dimension of B as declared
  98. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  99. *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
  100. *           be at least  max( 1, k ).
  101. *           Unchanged on exit.
  102. *
  103. *  BETA   - COMPLEX*16      .
  104. *           On entry, BETA specifies the scalar beta.
  105. *           Unchanged on exit.
  106. *
  107. *  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
  108. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  109. *           upper triangular part of the array C must contain the upper
  110. *           triangular part  of the  symmetric matrix  and the strictly
  111. *           lower triangular part of C is not referenced.  On exit, the
  112. *           upper triangular part of the array  C is overwritten by the
  113. *           upper triangular part of the updated matrix.
  114. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  115. *           lower triangular part of the array C must contain the lower
  116. *           triangular part  of the  symmetric matrix  and the strictly
  117. *           upper triangular part of C is not referenced.  On exit, the
  118. *           lower triangular part of the array  C is overwritten by the
  119. *           lower triangular part of the updated matrix.
  120. *
  121. *  LDC    - INTEGER.
  122. *           On entry, LDC specifies the first dimension of C as declared
  123. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  124. *           max( 1, n ).
  125. *           Unchanged on exit.
  126. *
  127. *
  128. *  Level 3 Blas routine.
  129. *
  130. *  -- Written on 8-February-1989.
  131. *     Jack Dongarra, Argonne National Laboratory.
  132. *     Iain Duff, AERE Harwell.
  133. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  134. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  135. *
  136. *
  137. *     .. External Functions ..
  138.       LOGICAL            LSAME
  139.       EXTERNAL           LSAME
  140. *     .. External Subroutines ..
  141.       EXTERNAL           XERBLA
  142. *     .. Intrinsic Functions ..
  143.       INTRINSIC          MAX
  144. *     .. Local Scalars ..
  145.       LOGICAL            UPPER
  146.       INTEGER            I, INFO, J, L, NROWA
  147.       COMPLEX*16         TEMP1, TEMP2
  148. *     .. Parameters ..
  149.       COMPLEX*16         ONE
  150.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  151.       COMPLEX*16         ZERO
  152.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  153. *     ..
  154. *     .. Executable Statements ..
  155. *
  156. *     Test the input parameters.
  157. *
  158.       IF( LSAME( TRANS, 'N' ) )THEN
  159.          NROWA = N
  160.       ELSE
  161.          NROWA = K
  162.       END IF
  163.       UPPER = LSAME( UPLO, 'U' )
  164. *
  165.       INFO = 0
  166.       IF(      ( .NOT.UPPER               ).AND.
  167.      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
  168.          INFO = 1
  169.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  170.      $         ( .NOT.LSAME( TRANS, 'T' ) )      )THEN
  171.          INFO = 2
  172.       ELSE IF( N  .LT.0               )THEN
  173.          INFO = 3
  174.       ELSE IF( K  .LT.0               )THEN
  175.          INFO = 4
  176.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  177.          INFO = 7
  178.       ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN
  179.          INFO = 9
  180.       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
  181.          INFO = 12
  182.       END IF
  183.       IF( INFO.NE.0 )THEN
  184.          CALL XERBLA( 'ZSYR2K', INFO )
  185.          RETURN
  186.       END IF
  187. *
  188. *     Quick return if possible.
  189. *
  190.       IF( ( N.EQ.0 ).OR.
  191.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  192.      $   RETURN
  193. *
  194. *     And when  alpha.eq.zero.
  195. *
  196.       IF( ALPHA.EQ.ZERO )THEN
  197.          IF( UPPER )THEN
  198.             IF( BETA.EQ.ZERO )THEN
  199.                DO 20, J = 1, N
  200.                   DO 10, I = 1, J
  201.                      C( I, J ) = ZERO
  202.    10             CONTINUE
  203.    20          CONTINUE
  204.             ELSE
  205.                DO 40, J = 1, N
  206.                   DO 30, I = 1, J
  207.                      C( I, J ) = BETA*C( I, J )
  208.    30             CONTINUE
  209.    40          CONTINUE
  210.             END IF
  211.          ELSE
  212.             IF( BETA.EQ.ZERO )THEN
  213.                DO 60, J = 1, N
  214.                   DO 50, I = J, N
  215.                      C( I, J ) = ZERO
  216.    50             CONTINUE
  217.    60          CONTINUE
  218.             ELSE
  219.                DO 80, J = 1, N
  220.                   DO 70, I = J, N
  221.                      C( I, J ) = BETA*C( I, J )
  222.    70             CONTINUE
  223.    80          CONTINUE
  224.             END IF
  225.          END IF
  226.          RETURN
  227.       END IF
  228. *
  229. *     Start the operations.
  230. *
  231.       IF( LSAME( TRANS, 'N' ) )THEN
  232. *
  233. *        Form  C := alpha*A*B' + alpha*B*A' + C.
  234. *
  235.          IF( UPPER )THEN
  236.             DO 130, J = 1, N
  237.                IF( BETA.EQ.ZERO )THEN
  238.                   DO 90, I = 1, J
  239.                      C( I, J ) = ZERO
  240.    90             CONTINUE
  241.                ELSE IF( BETA.NE.ONE )THEN
  242.                   DO 100, I = 1, J
  243.                      C( I, J ) = BETA*C( I, J )
  244.   100             CONTINUE
  245.                END IF
  246.                DO 120, L = 1, K
  247.                   IF( ( A( J, L ).NE.ZERO ).OR.
  248.      $                ( B( J, L ).NE.ZERO )     )THEN
  249.                      TEMP1 = ALPHA*B( J, L )
  250.                      TEMP2 = ALPHA*A( J, L )
  251.                      DO 110, I = 1, J
  252.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  253.      $                                          B( I, L )*TEMP2
  254.   110                CONTINUE
  255.                   END IF
  256.   120          CONTINUE
  257.   130       CONTINUE
  258.          ELSE
  259.             DO 180, J = 1, N
  260.                IF( BETA.EQ.ZERO )THEN
  261.                   DO 140, I = J, N
  262.                      C( I, J ) = ZERO
  263.   140             CONTINUE
  264.                ELSE IF( BETA.NE.ONE )THEN
  265.                   DO 150, I = J, N
  266.                      C( I, J ) = BETA*C( I, J )
  267.   150             CONTINUE
  268.                END IF
  269.                DO 170, L = 1, K
  270.                   IF( ( A( J, L ).NE.ZERO ).OR.
  271.      $                ( B( J, L ).NE.ZERO )     )THEN
  272.                      TEMP1 = ALPHA*B( J, L )
  273.                      TEMP2 = ALPHA*A( J, L )
  274.                      DO 160, I = J, N
  275.                         C( I, J ) = C( I, J ) + A( I, L )*TEMP1 +
  276.      $                                          B( I, L )*TEMP2
  277.   160                CONTINUE
  278.                   END IF
  279.   170          CONTINUE
  280.   180       CONTINUE
  281.          END IF
  282.       ELSE
  283. *
  284. *        Form  C := alpha*A'*B + alpha*B'*A + C.
  285. *
  286.          IF( UPPER )THEN
  287.             DO 210, J = 1, N
  288.                DO 200, I = 1, J
  289.                   TEMP1 = ZERO
  290.                   TEMP2 = ZERO
  291.                   DO 190, L = 1, K
  292.                      TEMP1 = TEMP1 + A( L, I )*B( L, J )
  293.                      TEMP2 = TEMP2 + B( L, I )*A( L, J )
  294.   190             CONTINUE
  295.                   IF( BETA.EQ.ZERO )THEN
  296.                      C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
  297.                   ELSE
  298.                      C( I, J ) = BETA *C( I, J ) +
  299.      $                           ALPHA*TEMP1 + ALPHA*TEMP2
  300.                   END IF
  301.   200          CONTINUE
  302.   210       CONTINUE
  303.          ELSE
  304.             DO 240, J = 1, N
  305.                DO 230, I = J, N
  306.                   TEMP1 = ZERO
  307.                   TEMP2 = ZERO
  308.                   DO 220, L = 1, K
  309.                      TEMP1 = TEMP1 + A( L, I )*B( L, J )
  310.                      TEMP2 = TEMP2 + B( L, I )*A( L, J )
  311.   220             CONTINUE
  312.                   IF( BETA.EQ.ZERO )THEN
  313.                      C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
  314.                   ELSE
  315.                      C( I, J ) = BETA *C( I, J ) +
  316.      $                           ALPHA*TEMP1 + ALPHA*TEMP2
  317.                   END IF
  318.   230          CONTINUE
  319.   240       CONTINUE
  320.          END IF
  321.       END IF
  322. *
  323.       RETURN
  324. *
  325. *     End of ZSYR2K.
  326. *
  327.       END
  328.